// Gmsh project created on Tue Apr 11 15:50:36 2017
// Merge "dbWide.brep";
// To generate mesh:\n
// gmsh -3 -optimize_ho -format msh -o dbWide.msh dbWide.geo

// Units are dimensionless, but need to be
// consistent. Consider millimeters

// Convert mm to metres /////////////////////////
mm2m = 0.001;
// Point tolerance //////////////////////////////
lc = 10e-3 * mm2m;
// Origin ///////////////////////////////////////
orgn_x = 0;
orgn_y = 0;
orgn_z = 0;
// Thickness ////////////////////////////////////
thick = 0.5 * mm2m;
// // Total height /////////////////////////////////
// total_hei = 12.5 * mm2m;

/////////////////////////////////////////////////
//                Mid substance                //
/////////////////////////////////////////////////
// ** Mid substance width
mid_subst_wid = 0.6 * mm2m;
// ** Mid substance height
mid_subst_hei = 6 * mm2m;
// ** Mid substance number of points in x
mid_subst_x_div = 7;
// ** Mid substance number of points in z
mid_subst_z_div = mid_subst_x_div * 3;
// ** Mid substance bottom line
// *** Points of the bottom line of the mid
// substance
mid_subst_P1 = newp;
Point(mid_subst_P1) = {orgn_x, orgn_y, orgn_z, lc};
mid_subst_P2 = newp;
Point(mid_subst_P2) = {mid_subst_wid, orgn_y, orgn_z, lc};
// *** Line
mid_subst_L1 = newl;
Line(mid_subst_L1) = {mid_subst_P1, mid_subst_P2};
// ** Mid substance top points
Translate {0, 0, mid_subst_hei} {
  Duplicata { Point{mid_subst_P1, mid_subst_P2}; }
}
/* The new two points from the Duplicata have the
index of mid_subst_P2 +1 and +2 (mid_subst_P2 is
the last point so far) */
mid_subst_P3 = mid_subst_P2 + 1;
mid_subst_P4 = mid_subst_P3 + 1;
// ** Mid substance vertical and top lines
mid_subst_L2 = newl;
Line(mid_subst_L2) = {mid_subst_P1, mid_subst_P3};
mid_subst_L3 = newl;
Line(mid_subst_L3) = {mid_subst_P3, mid_subst_P4};
mid_subst_L4 = newl;
Line(mid_subst_L4) = {mid_subst_P4, mid_subst_P2};
// ** Join lines for mid substance
mid_subst_Lp1 = newll;
Line Loop(mid_subst_Lp1) = {mid_subst_L1,
  -mid_subst_L4, -mid_subst_L3, -mid_subst_L2};
// *** Print line number (example)
Printf("Line number mid_subst_Lp1: ’%g’",
  mid_subst_Lp1);
// ** Discretise lines (for a controlled mesh)
Transfinite Line{mid_subst_L1, mid_subst_L3} = mid_subst_x_div;
Transfinite Line{mid_subst_L2, mid_subst_L4} = mid_subst_z_div Using Bump 0.3;
// ** Create surface
mid_subst_S1 = news;
Plane Surface(mid_subst_S1) = {mid_subst_Lp1};
// *** Discretise surface by specifying the four
//corners of the interpolation
// forces meshing to contain structured triangles
Transfinite Surface{mid_subst_S1} = {mid_subst_P1,
  mid_subst_P3, mid_subst_P4, mid_subst_P2};
// Hexas instead of tetras
Recombine Surface{mid_subst_S1};

/////////////////////////////////////////////////
//                   Chamfer                   //
/////////////////////////////////////////////////
// * Bottom line of chamfer
// ** Reuse top points of mid-substance as bottom
//    points of the chamfer
chmf_P1 = mid_subst_P4;
chmf_P2 = mid_subst_P3;
// ** Reuse top line of mid-substance as bottom
//    line of chamfer
chmf_L1 = mid_subst_L3;
// * Arc of chamfer
// ** Radius
chmf_R = 3.5 * mm2m;
// ** Centre of arc (same height as bottom line;
//    it needs to be tangent to the vertical line
//    of the mid-substance).
chmf_Cx = chmf_R + mid_subst_wid;
chmf_Cz = mid_subst_hei;
chmf_C = newp;
Point(chmf_C) = {chmf_Cx, orgn_y, chmf_Cz, lc};
// ** End point of arc
// *** To get the position of this point, go to
//     FreeCad, find the index # of the point
//     (in the sketch) and do:
// App.ActiveDocument.Sketch014.getPoint(#, 1)
//     I'm sure that you can also calculate
//     the tangency conditions of the two arcs
//     comprising the chamfer and clamping
//     regions ;) ...
chmf_P3_x = 2.323068841391659 * mm2m;
chmf_P3_z = 9.015359609265856 * mm2m;
chmf_P3 = newp;
Point(chmf_P3) = {chmf_P3_x, orgn_y, chmf_P3_z};
// ** Arc (line) of chamfer
chmf_A = newl;
Circle(chmf_A) = {chmf_P1, chmf_C, chmf_P3};
// * Top line of chamfer
// ** End point
chmf_P4 = newp;
Point(chmf_P4) = {orgn_x, orgn_y, chmf_P3_z};
// ** Top line of chamfer
chmf_L3 = newl;
Line(chmf_L3) = {chmf_P3, chmf_P4};
// * Vertical line of chamfer
chmf_L4 = newl;
Line(chmf_L4) = {chmf_P4, chmf_P2};
// * Closed line for chamfer
chmf_Lp1 = newll;
Line Loop(chmf_Lp1) = {chmf_L1, chmf_A,
  chmf_L3, chmf_L4};

// ** Discretise lines (for a controlled mesh)
// *** Number of points in x
chmf_x_div = mid_subst_x_div;
chmf_z_div = mid_subst_x_div * 2;
Transfinite Line{chmf_L1, chmf_L3} = chmf_x_div;
Transfinite Line{chmf_L4} = chmf_z_div Using Bump 0.8;
Transfinite Line{chmf_A} = chmf_z_div Using Progression 1.046;
// ** Create surface
chmf_S1 = news;
Plane Surface(chmf_S1) = {chmf_Lp1};
// *** Discretise surface by specifying the four
//corners of the interpolation
// forces meshing to contain structured triangles
// Transfinite Surface{chmf_S1};  // works
Transfinite Surface{chmf_S1} = {chmf_P1,
  chmf_P2, chmf_P4, chmf_P3};
// Hexas instead of tetras
Recombine Surface{chmf_S1};

/////////////////////////////////////////////////
//                   Extrude                   //
/////////////////////////////////////////////////
Extrude {0, thick, 0} {
  // Use the mid-substance, and the chamfer (no clamp)
  Surface{mid_subst_S1, chmf_S1};
  // Four layers
  Layers{4};
  // I wished that PreView could do hexas with Gmsh
  // Comment this block if using FEBio
  Recombine;
}
// With this line, PreView can import the mesh
Physical Volume("dbWide_vol") = {1, 2};

// Comment these lines when exporting to FEBio PreView
// (got numbers with help of the GUI: Tools -> Visibility ->
// Tree browser -> click on a surface -> Apply)
Physical Surface("bottom") = {20};
Physical Surface("back") = {mid_subst_S1, chmf_S1};
Physical Surface("left") = {32, 54};
Physical Surface("top") = {50};
Physical Surface("front") = {33, 55};
